getwd()

rm(list = ls())
# install.packages('ggplot2')

library(mgcv) 
library(ggplot2)
library(reshape2) 
library(tidyverse)
source("toolbox.R")


# parameter
mech <- 5 # choose 1,2,3,4 for other mechanism
N_sim <- 500
alpha <- 0
beta <- 0
a <- 0.31
b <- 0.69

# 描述不同fitness定义方式
# for antagonist
F_antagonist <- function(mech){
  E1 <- (1-H_V_A) # n=1 for A_min & P_max
  E2 <- (1-H_V_A)-alpha*(1-H_V_Po)
  E3 <- (1-H_V_A)
  E4 <- (1-H_V_A)+beta*(1-H_V_Po)
  E5 <- (1-H_V_A)+beta*(1-H_V_Po)
  
  E <- list(E1,E2,E3,E4,E5)
  return(E[[mech]])
}

# for pollinators
F_pollinators <- function(mech){
  E1 <- (1-H_V_Po)   # n=1 for A_min & P_max
  E2 <- (1-H_V_Po)-alpha*(1-H_V_A)
  E3 <- (1-H_V_Po)-alpha*(1-H_V_A)
  E4 <- (1-H_V_Po)-alpha*(1-H_V_A)
  E5 <- (1-H_V_Po)+alpha*(1-H_V_A)
  
  E <- list(E1,E2,E3,E4,E5)
  return(E[[mech]])
}

# for plant
F_plant <- function(mech){
  E1 <- (a*H_A_V+b*(1-H_Po_V))
  E2 <- (a*H_A_V+b*(1-H_Po_V))
  E3 <- (a*H_A_V+b*(1-H_Po_V))
  E4 <- (a*H_A_V+b*(1-H_Po_V))
  E5 <- (a*H_A_V+b*(1-H_Po_V))
  
  E <- list(E1,E2,E3,E4,E5)
  return(E[[mech]])
}



AP_obs  <- read.csv("data/AP_Africa_B.csv", header = TRUE, as.is = TRUE, row.names = 1)
PoP_obs <- read.csv("data/PoP_Africa.csv", header = TRUE, as.is = TRUE, row.names = 1)
PV_obs  <- read.csv("data/PV_Africa.csv", header = TRUE, as.is = TRUE, row.names = 1)

index_choose <- c(35,40,16,40,74,12,14,15,1,19) # 按照原则一抽取的 n 种 VOCs 
PV_obs_1 <- PV_obs[,c(index_choose)]
index_all <- 1:ncol(PV_obs)

PV_obs_2 <- PV_obs[,c(index_all[-index_choose])] 
PV_obs_3 <- PV_obs_2[,c(colSums(PV_obs[,index_all[-index_choose]])<=3)]

index_1 <- 1:ncol(PV_obs_3)
# index_ <- sample(index_1 , 56)
# print(index_)
# PV_obs_4 <- PV_obs_3[,c(index_ )]

PV_obs <- cbind(PV_obs_1,PV_obs_3)



nA  <- nrow(AP_obs) # 植食动物的数量 （行计数）
nP  <- ncol(AP_obs) # 植物的数量 （列计数）
nV  <- ncol(PV_obs) #  VOC 的数量
nPo <- nrow(PoP_obs) # 传粉蜂的数量


# index_ <- sample(index , 22)


# x <- combn(index,30)

# PV_obs<-PV_obs[,c(index_ )] #取数据集第1到3列
## 生成模拟的初始化矩阵
# 1、初始矩阵不影响模拟结果
# 2、初始时要保证任意植物至少有与一种植食动物有关；任意VOC都至少由一种植物释放。因此下述的colSum 的所以值都大于 0
AP  <- matrix(rbinom(nA*nP,1,0.5),nA,nP)
PV  <- matrix(rbinom(nP*nV,1,0.5),nP,nV)
PoP <- matrix(rbinom(nPo*nP,1,0.5),nPo,nP)
colSums(AP) #  确保 colSum 的所以值都大于 0
colSums(PV) 
colSums(PoP)
# 每一次循环基因变异的比例，该值仅影响收敛速度，不影响结果。
M1 <- 0.2*nP*nV/2 # PV 矩阵突变的比例
M2 <- 0.2*nA*nP # AP 矩阵突变的比例
M3 <- 0.2*nPo*nP # PoP 矩阵突变的比例

#######################################
### Analyses1: Simulation process #####
#######################################
# create plot
par(mar=c(4,4,2,9))
plot(-1, xlim = c(0,N_sim), ylim = c(0,1), ylab = "Fitness", xlab = "Time")
legend(x=N_sim+N_sim/50, y=0.8, title = "Simulated", legend = c("Plant Fitness","Antagonist Fitness","pollinators Fitness"),
       pch = c(3,4,5), col = c("green3","red","blue"), xpd=TRUE, bty="n",title.font = 2) 


# 根据实地观察，计算并绘制熵值和fitness
{
  AV_obs  <- as.matrix(AP_obs) %*% as.matrix(PV_obs)
  PoV_obs <- as.matrix(PoP_obs) %*% as.matrix(PV_obs)
  # "H_A_VOC" 函数 i:矩阵 o:互信息/条件熵，toolbox
  # H_V   <- H_A_VOC(PV_obs)[["Hn_V"]]
  # H_P_V <- H_A_VOC(PV_obs)[["Hn_S_V"]]
  # H_V_P <- H_A_VOC(PV_obs)[["Hn_V_S"]]
  H_A_V <- H_A_VOC(AV_obs)[["Hn_S_V"]]
  H_V_A <- H_A_VOC(AV_obs)[["Hn_V_S"]]
  # H_A_P <- H_A_VOC(AP_obs)[["Hn_S_V"]]
  # H_P_A <- H_A_VOC(AP_obs)[["Hn_V_S"]]
  H_Po_V <- H_A_VOC(PoV_obs)[["Hn_S_V"]]
  H_V_Po <- H_A_VOC(PoV_obs)[["Hn_V_S"]]
  # H_Po_P <- H_A_VOC(PoP_obs)[["Hn_S_V"]]
  # H_P_Po <- H_A_VOC(PoP_obs)[["Hn_V_S"]]
  
  # 对抗系统中的 fitness 表达, toolbox
  E_plant       <- F_plant(mech) 
  E_antagonist  <- F_antagonist(mech) 
  E_pollinators <- F_pollinators(mech) 
  
  
  
  # E_obs <- c(0, E_plant, E_antagonist, E_pollinators, H_V, 
  #            H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
  # names(E_obs) <- c("N","E_plant", "E_antagonist", "E_pollinators", "H_V", "H_P_V", "H_V_P", "H_A_V", "H_V_A", "H_A_P", "H_P_A","H_Po_P","H_P_Po","H_Po_V","H_V_Po")
  # 
  ## 野外实际观测值划线，图A
  # 条件熵 for PV, AV and AP
  # abline(h = H_P_V, lty = "dotted", pch = 8)
  # abline(h = H_A_V, lty = "solid", pch = 8)
  # abline(h = H_A_P, lty = "dashed", pch = 8) 
  ## 植物和植食动物的fitness，图B
  abline(h = E_antagonist, col = "red")
  abline(h = E_pollinators, col = "blue")
  abline(h = E_plant, col = "green3") # the same as H(A|V) in mech = 1
}
